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Finite  Methods  for  a  Nonlinear  Allocation  Problem 

Alan  Washburn 


1 .  Introduction 

We  consider  here  finite  algorithms  for  solving  the  nonlinear  program  PI: 

minimize  f(yi,...,yn) 

m 

subject  to  X 

i=l 
n 

X  Xij=bi;  l^<m, 

j=i 

xij>0;  l<i^,  l<j^. 


f\-\ 


It  is  assumed  in  PI  that  bi^  for  all  i  and  that  E(ij)^  for  all  (i ,j),  with  at  least  one  positive 
E(i  j)  in  each  row  i.  It  is  further  assumed  that  f(y)  is  continuous  and  decreasing  in  each  of 
its  arguments  on  some  convex  set  S  that  includes  the  nonnegative  orthant  of  R".  y 
represents  the  vector  (yi,...,yn)i  sJ'd  similarly  x  will  represent  the  collection  of  all  mn  of 
the  Xjj,  etc.  PI  is  a  special  case  of  the  nonlinear  network  problem,  a  problem  that  is  already 
known  to  possess  sufficient  special  structure  to  make  the  development  of  special  purpose 
software  attractive(Ahlfeld,  Dembo,  Mulvey,  and  Zenios(1987)).  Our  intention  is  to 
specialize  PI  even  further  by  making  restrictive  assumptions  about  f(y).  The  components 
of  y  will  be  called  "potentials",  one  for  each  column  of  the  matrix  (E(i,j)).  Potentials  are 
not  logically  necessary  in  defining  PI,  since  the  expressions  defining  y  could  simply  be 
substituted  into  f(y),  but  the  generous  notation  will  prove  useful  in  the  sequel.  The 
components  of  x  will  be  call^  "allocations." 

PI  could  be  interpreted  as  a  Search  Theory  problem  by  assuming  that  search  is  to  be 
conducted  over  a  fixed  time  period  by  m  distinct  types  of  searcher,  with  bj  units  of  type  i 
effort  available  over  the  period.  Assume  that  a  single  target  is  located  with  probability  pj  in 
one  of  n  regions,  that  region  j  has  area  Aj,  and  that  the  reward  for  finding  the  target  in 
region  j  is  Vj.  Assume  further  that  the  amount  of  area  swept  by  a  searcher  of  type  i  in 
region  j  is  sjj.  If  Xjj  is  the  allocation  of  search  effort  of  type  i  to  region  j,  then  the  total  area 

m 

swept  in  region  j  is  X  xijSij.  If  the  search  in  each  area  is  "random"  (Koopman(1956)), 

i=l 

then  the  number  of  times  a  target  in  region  j  is  detected  is  a  Poisson  random  variable  with 


=yj  (so  let  E(i,j)=Sij/Aj),  and  the  probability  of  no  detections  is 


n 

exp(-yj).  The  average  "regret"  in  terms  of  value  not  found  is  then  f(y)=X  ^jPj®’‘P(~yj)- 

j=i 

PI  is  then  the  problem  of  allocating  m  resources  to  n  regions  in  order  to  minimize  the 
average  undetected  value.  PI  can  be  thought  of  as  a  generalization  from  the  case  m=l, 
which  already  has  specially  tailored  solution  procedures  (Washbum(1981)). 

The  US  Air  Force's  Heavy  Attack  model  (Clasen,  Graves,  and  Lu(1974))  is  also 
closely  related  to  PI ,  with  resources  being  aircraft  sonies,  areas  being  classes  of  targets. 
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and  potentials  being  "average  number  of  targets  killed  were  it  not  for  the  fact  that  some  of 
the  potential  will  be  'wasted'  in  attacking  targets  already  dead".  In  general,  PI  is  applicable 
when  n  projects  must  be  simultaneously  undertaken  by  employing  m  resources,  y  is  a 
vector  measure  of  progress,  and  f(y)  is  some  global,  scalar  measure  of  work  not 
completed. 

2 .  Basic  Feasible  Solutions 

PI  has  mn  potentially  positive  allocations,  but  optimal  solutions  invariably  have  most 
of  these  being  0.  The  reason  for  this  is  that  optimal  solutions  need  never  include  cycles  of 
positive  allocations. 

Definition:  If  a  feasible  allocation  for  PI  has  an  alternating  sequence 
jbil02.i2.-”.iLjL+l  of  columns  and  rows  all  distinct  except  that  jL+i=ji,  and  if  Xij>0  for 
every  consecutive  pair  in  the  sequence  ((i,j)=(ikJk)  or  (i,j)=(ik,  jk+i)  for  some  k^),  then 
the  sequence  is  a  cycle. 

Definition:  A  feasible  allocation  for  PI  is  conservative  if  xij=0  whenever  E(i,j)=0. 

Theorem  1:  Given  the  assumptions  of  section  1,  PI  has  an  optimal,  conservative 
solution  with  no  cycles. 

Proof:  PI  has  an  optimal  solution  because  the  objective  function  regarded  as  a  function 
of  X  is  continuous  on  a  compact  set  (Bazaraa  and  Shetty(1979)).  An  optimal  solution  x 
can  easily  be  converted  to  a  conservative  optimal  solution  by  shifting  any  offending 
allocation  xjj  from  column  j  to  some  other  column  k  where  E(i,k)>0  (by  assumption  at 
least  one  sucti  column  k  always  exists  for  every  i).  Assume  therefore  that  x  is  optimal  and 
conservative,  but  that  a  cycle  exists.  In  that  case  construct  a  modified  solution  x'  by  adding 

5(i,j)  to  Xjj  for  each  consecutive  pair  (i,j)  in  the  sequence,  where  6(i,j)  is  defined  as 
follows: 

5(ii,ji)=z  (an  arbitrary  real  number) 

5(ik0k+l)=- 6(ik.jk);  l^^L 

5(ik+l0k+l)=-  5(ikOk+l)E(ikOk+l)/E(ik+lOk+l);  l^^-l. 

n 

The  numbers  5(i,j)  are  constructed  so  that  ^  x'ij=bi  for  all  i  as  before.  Also 

j=i 

m 

^  x'ijE(i,j)=yj  as  before,  except  that  the  latter  sum  for  j=ji  is  increased  by 
i=i 

5(ii,ji)E(ii,ji)-5(iL,ji)E(iL,ji).  This  difference  is  proportional  to  z.  If  the  proportionality 
constant  is  nonnegative,  increase  z  until  x'ij=0  for  some  (i,j)  in  the  cycle;  there  will  be 

some  such  (i,j)  because  the  numbers  5ij  are  proportional  to  z  with  nonzero  proportionality 
constants  that  alternate  in  sign.  Otherwise,  decrease  z  to  achieve  the  same  end.  Since  f()  is 
decreasing  in  each  variable,  the  x'  solution  is  at  least  as  good  as  the  x  solution,  with  one 
less  positive  variable.  By  repeating  this  operation  if  necessary,  an  optimal  solution  with  no 
cycles  can  be  found.  QED 

Since  an  optimal  conservative  solution  with  no  cycles  exists,  it  makes  sense  to  search 
for  the  optimal  solution  amongst  solutions  with  that  property.  A  major  advantage  of  such 
solutions  is  that  there  can  be  at  most  n+m-1  positive  variables  Xjj.  To  see  this,  consider  the 
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nonoriented  bipartite  graph  G  formed  by  connecting  m  "row  nodes"  to  n  "column  nodes", 
with  i  being  connected  to  j  if  and  only  if  Xij>0.  Such  a  graph  must  consist  of  a  number  of 
connected  components  Tj,. .  .Tk,  each  of  which  is  not  connected  to  any  of  the  others. 
Some  of  these  components  may  contain  a  single  node  and  no  edges,  but  every  node  is  in 
exactly  one  component.  Since  each  component  is  connected  and  has  no  cycles,  it  is  a  tree 
(it  is  convenient  to  refer  to  even  singleton  components  as  "trees"),  so  G  is  a  "forest".  If 
has  He  nodes,  it  must  have  exactly  r^-l  edges  (Berge(1962)),  even  if  rk=l.  Therefore  the 
total  number  of  edges  is  n+m-K,  which  is  at  most  n+m-1. 

The  graph  G  described  above  will  be  called  the  "basis"  of  any  conservative  solution 
with  no  cycles.  More  generally. 

Definition:  Any  nonoriented  bipartite  graph  G  with  the  following  properties  is  a  basis: 

1)  The  nodes  of  G  are  all  of  the  m+n  row  and  column  nodes. 

2)  All  row  nodes  I  for  which  bi>0  have  at  least  one  incident  edge. 

3)  If  E(i,j)=0,  then  there  is  no  edge  (i,j). 

4)  There  are  no  cycles. 


The  requirement  that  a  solution  should  have  no  cycles  is  sufficient  to  establish  that  the 
corresponding  basis  is  a  forest,  but  it  is  not  sufficient  to  establish  a  one-to-one  relationship 
between  bases  and  solutions.  In  order  to  establish  such  a  relationship,  further 
assumptions  about  f()  are  needed. 

Definition:  f()  is  good  if  it  has  the  following  properties: 

1)  f()  is  differentiable  and  decreasing  on  the  interior  of  the  convex  set  S  described 
earlier.  Let  f()  be  the  gradient  of  f(). 

2)  There  exists  a  unique  function  g(iJ.)  taking  values  in  S  such  that  f(g(|i))-(-|i=0  for 
all  p>0,  where  by  p>0  we  mean  |ij>0  for  j=l,...,n. 

3)  If  C  is  a  nonempty  subset  of  { l,...,n),  if  ^>0,  and  if  scalar  p>0,  then  the  equation 

]^lijgj(z|i)=p  has  a  unique  positive  solution  z.  Furthermore  z  is  a  continuous  function  of 
jec 


p  for  p>0. 

In  practice  the  methods  introduced  below  will  be  attractive  only  if  the  solution  z 
required  in  the  last  part  is  easily  computed,  since  the  equation  will  have  to  be  solved  many 
times.  Here  are  two  examples  of  good  functions: 

n 

Example  1:  f(y)=X  S=R".  This  is  the  Search  Theory  function 

j=i 

introduced  earlier,  except  that  pj  has  been  incorporated  into  Vj.  In  this  case  gj(|i.)=ln(Vj/pj) 


for  j=l,...,n,  and  ln(z)= 


Xpjgj(p)-p 

Uec 


LjeC  J 


Example  2:  f(y)=X  l/(l+yj),  with  S={y:  yj>-l  for  j=l,...,n}.  Here  gj()i)=(Vvj/|ij  -l) 
j=i 


and 


J£C 


JEC  . 


,  well  defined  for  p>0. 


such  that 


The  Kuhn-Tucker  conditions  for  PI  are  that  there  should  exist  )i,|i,x,  and  y 
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KTl)  A,i>|ijE(i  j)  for  all  ij 
KT2)  xij>Oforalli,j 

KT3)  X,i=^.jE(iJ)  when  Xij?tO 

m 

KT4)  Y,  ^ijE(io)=yj  all  j 

i=i 

n 

KT5)  ^  xij=bi  for  all  i,  and 
j=i 

KT6)  y=g(^). 


The  object  now  is  to  use  the  equality  parts  of  the  KT  conditions  (KT3-KT6)  to  set  up  a 
one-to-one  correspondence  between  bases  and  solutions  for  good  objective  functions,  after 

which  solutions  corresponding  to  bases  will  be  called  "basic."  The  vectors  X  and  |i  will 
each  be  called  "multipliers."  Basic  solutions  will  have  xij=0  for  nonbasic  edges  (i,j),  but 
the  following  theorem  permits  one  exception  to  support  the  pivoting  ideas  of  section  4. 


Theorem  2:  For  any  tree  T  let  R(T)  and  C(T)  be  the  row  nodes  and  column  nodes  of  T, 


respectively.  Let  G  be  a  basis  composed  of  trees  (Ti,...,Tk),  and  suppose  that  Xi=pjE(i,j) 
for  all  edges  (i,j)  of  G.  Then 

a)  If  (x,y)  solve  KT4-KT5,  and  if  all  nonbasic  allocations  except  for  xjj  are  zero,  then 
for  any  component  T  of  G, 


S  iijyj 

JeC(T) 


-WxuE(I,J)5(C(T),J) 


ieR(T) 


-Xixu5(R(T),I), 


(1) 


where  now  6(W,w)  is  1  if  w  is  in  the  set  W,  otherwise  0.  Furthermore, 

b)  If  (X,,p.,y)  is  such  that  (1)  holds  for  every  component  T  of  G,  then  KT4-KT5  have 
a  unique  solution  x  for  which  all  nonbasic  allocations  except  xn  are  zero. 


Proof  of  part  a)  Multiply  both  sides  of  KT4  by  pj  and  sum  to  obtain 

m 

Z  Z  l^j^ijE(i.j)=  X  l^jyj- 

jeC(T)  >=1  j£C(T) 

Since  xij=0  when  j£C(T)  unless  either  i=I  or  ieR(T), 

Z  X  ^jXijE(i,j)  +  pjxijE(l,J)5(c(T),J)  =  pjyj. 

jEC(T)  i£R(T)  jEC(T) 

But  |ijE(i,j)=?ii  when  jeC(T)  and  ieR(T)  by  assumption,  so 

Z  Z’^iJ  +  ltJxijE(I,J)5(c(T),J)  =  Zl^jyj- 

i£R(T)  j£C(T)  j£C(T) 

Since  KT5  holds, 

Z  =bi  -  xu5({i},D  for  all  ieR(T). 

jEC(T) 

Part  a)  follows  upon  substituting  (5)  into  (4). 


(2) 


(3) 

(4) 

(5) 
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Proof  of  part  b):  To  avoid  unnecessary  complication  assume  xu=0;  otherwise  bj  and  yj 
can  be  adjusted  by  subtracting  xu  and  xuE(IJ),  respectively.  Consider  any  component  T 
of  G.  We  must  show  that  there  is  exactly  one  way  to  assign  allocations  Xij  to  the  edges  of 
T  that  is  consistent  with  KT4-KT5.  This  is  trivial  if  T  is  a  singleton;  otherwise,  since  T  is 
a  tree,  there  must  be  at  least  1  pendant  (connected  to  exactly  one  other  node)  node.  If  this 
pendant  node  is  row  P,  let  row  P  be  connected  to  column  Q.  Then  xpQ  is  the  only  nonzero 
allocation  in  row  P,  and  is  therefore  necessarily  bp  by  KT5.  Define  y'  by 
y'Q=yQ-  bpE(P,Q),  y  j=yj  for  J^tQ,  and  let  T'  be  the  tree  remaining  after  node  P  and  edge 

(P,Q)  are  deleted.  Since  Xp=|iqE(P,Q),  |iQy'Q=|iQyQ-Xpbp.  Therefore 

X  J  “  X  pendant  node  is  column  Q,  let  it  be  connected  to  row  P. 

jeC(T')  icRO  ■; 

Since  xpq  is  the  only  positive  allocation  in  column  Q,  necessarily  xpq  is  yQ/E(P,Ci) 
according  to  KT4.  xpg  is  well  defined  because  E(P,Q)>0.  Let  b'p=bp-xpQ,  otherwise  let 
b'i=bi  for  i?tP,  and  let  T'  be  the  tree  resulting  from  deleting  column  Q  and  edge  (P,Q). 

Since  Xp=)IqE(P,Q),  Xpb'p=Xpbp-|iQyQ.  Therefore  ^  |j,jyj  =  ^  ^ib';.  Since  in 

j£C(T)  i£R(T') 

either  case  T'  has  one  less  node  than  T,  and  since  the  proposition  is  true  for  graphs  with 
one  node,  the  theorem  is  proved  by  induction.  QED 

Multipliers  for  which  Xi=|ijE(i,j)  on  basic  edges  are  easy  to  determine  if  there  is  no 
requirement  for  equation  (1)  to  hold.  The  reason  is  that  if  A.;  or  |ij  is  determined  at  any 

node,  then  the  same  thing  is  true  of  all  neighbors  of  the  node.  Thus  X  and  |i  can  be 
detennined  in  any  component  T  of  G  by  assigning  an  arbitrary  positive  value  to  any  node. 
Let  (A*,|i*)  be  such  multipliers.  Then  any  other  set  (A,,p)  with  the  same  propeny  is  some 
scalar  multiple  of  (A*,p*).  Let  z>0  be  the  multiple.  Then  z  must  be  determined  so  that  (1) 
holds.  After  substituting  gj(z^*)  for  yj,  (1)  is: 

X  -li*jxijt.(l,J)5(C(T),J)=  X  -  ?i*ixij6(R(T),l).  (6) 

jeC(T)  ieRCT) 

As  long  as  xij<bi,  the  right  hand  side  of  (6)  is  nonnegative.  Equation  (6)  will  then  have  a 
unique  positive  solution  z  as  long  as  f()  is  good.  Given  z,  KT^KT6  have  the  unique 

solution  (A.,|i.)=z(A.*,|i*),  y=g(M.),  and  x  the  unique  solution  of  KT4-KT5.  The  x  part  is 
the  desired  basic  solution  corresponding  to  the  basis  G-hereafter  the  G-basic  solution. 

Thus  the  unique  G-basic  solution  can  be  easily  determined  as  long  as  xu<bi . 

A  G-basic  solution  need  not  satisfy  KT2,  but  if  it  does  then  it  will  be  called  a  basic 
feasible  solution  and  G  will  be  called  "feasible".  Some  basic  feasible  solution  is  optimal, 
so  PI  could  be  solved  in  principle  by  examining  all  bases.  In  analogy  with  the  Simplex 
Method,  however,  the  object  is  to  avoid  this  by  utilizing  a  procedure  that  is  directional  in 
the  sense  of  considering  a  sequence  of  feasible  bases,  each  of  which  is  better  than  the  last. 
Two  of  these  will  be  considered  below.  In  each  case  we  will  prove  only  that  the  objective 
function  sequence  is  nonincreasing,  rather  than  decreasing,  with  equality  being  possible  in 
degenerate  cases.  Rather  than  deal  with  the  distracting  issue  of  degeneracy  here,  we  simply 
assume  nondegeneracy  for  the  next  two  sections.  There  is  every  reason  to  expect  that  the 
same  techniques  that  handle  degeneracy  in  the  Simplex  Method  will  also  work  here. 
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3 .  The  Manifold  Suboptimization  Method 

In  this  section  nonbasic  variables  will  always  be  0. 

A  G-basic  solution  is  the  optimal  solution  of  P2(G): 

mimimize  f(y) 

m 

subject  to  E(i,j)xij=yj;  l<j<n, 

i=l 
n 

^  xij<bi;  l^<m, 
j=i 

xij=0  unless  (i  j)eG. 

To  prove  this  note  that  the  G-basic  solution  solves  KT3-KT6.  But  KT3-KT6  are  the 
Kulm  Tucker  conditions  for  P2(G),  and  the  Kuhn  Tucker  conditions  are  sufficient  for  a 
global  minimum  because  f()  is  convex  as  a  function  of  x.  P2(G)  does  not  require 
allocations  to  be  nonnegative,  so  the  feasible  space  is  not  compact.  Absent  the  constraints 
forcing  nonbasic  variables  to  be  0,  P2(G)  might  lack  a  (finite)  optimal  solution. 
Nonetheless,  the  G-basic  solution  is  optimal. 

Given  a  basis  G  and  a  corresponding  solution  x  feasible  in  PI,  suppose  KTl  is  not 
satisfied  for  some  (I,J).  The  natural  computational  step  is  to  simply  put  (I,J)  into  G.  Any 
cycle  that  forms  (a  cycle  will  form  if  and  only  if  I  and  J  are  in  the  same  tree)  can  always  be 
profitably  destroyed  by  removing  some  other  edge  (Ki,Li)  determined  as  in  Theorem  1. 
(I,J)  would  be  (ii,ji)  in  that  theorem,  and  z  would  be  positive  because  the  reduced  cost  for 
(I,J)  is  negative.  Let  G^  be  G  with  (I,J)  added  and,  if  necessary,  (Kid^j)  deleted.  Let  x^+ 
be  X  if  (Kid-i)  is  not  deleted,  or  otherwise  the  solution  x'  of  Theorem  1.  xl+  is  better  than 
X  and  feasible,  but  is  not  in  general  the  G*-basic  solution.  The  latter  solution  (call  it  x^)  is 
better  than  x^+  because  it  is  optimal  in  P2(G*)  while  x^+  is  merely  feasible,  x^  is  therefore 
basic  and  better  than  x.  However,  x^  may  not  be  feasible  in  PI. 

If  x^  has  one  or  more  negative  components,  let  x2+=x(l-a)+xla,  where  of  [0,1]  is 
selected  so  that  x2+  is  the  closest  nonnegative  point  to  xL  Let  (K2,L2)  be  an  edge  in  G^ 
whose  allocation  is  driven  to  0  at  x2+.  Since  f()  is  convex  and  is  better  than  x,  x2+  is 
also  better  than  v  x2+  is  also  feasible  in  P2(G2),  where  G^  ic  the  basis  obtained  by 
deleting  edge  (K2,L2)  from  G^  Let  x^  be  optimal  in  G^,  so  that  x^  is  better  than  x2+.  If 
x2  also  has  negative  components,  continue  to  define  x^+.G^,  x^,  x^'*',  etc.,  until  finally  x*^ 
is  nonnegative,  as  must  eventually  happen  because  an  edge  is  deleted  from  the  basis  at  each 
step,  x*'  is  then  basic,  feasible  (in  PI),  and  (barring  degeneracy)  better  than  x.  Since  there 
are  only  finitely  many  bases,  none  of  which  can  repeat  because  the  objective  function 
decreases  at  each  step,  the  optimum  solution  must  eventually  be  found  in  firdtcly  many 
steps. 

The  above  algorithm  is  sufficiently  similar  to  Zangwill's  (1970)  method  of  Manifold 
suboptimization  that  the  name  has  been  retained,  but  there  is  a  difference.  The  difference  is 
that  in  Zangwill's  method  a  G-basic,  feasible  solution  that  is  not  optimal  would  lead  to  next 
considering  a  modified  basis  G'=G+(I,J),  where  (I,J)  is  some  edge  at  which  the  reduced 
cost  is  negative.  Unfortunately  G'  may  have  a  cycle,  in  which  case  P2(G')  may  not  have  a 
finite  optimal  solution.  The  method  described  above  works  because  every  manifold 
considered  has  a  basis  that  is  a  forest  of  trees. 

Manifold  suboptimization  essentially  operates  by  increasing  xjj  so  much  that  the  next 
basic  solution  may  not  be  feasible,  and  then  fixing  the  problem.  If  only  feasible  solutions 
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are  to  be  considered,  then  a  test  for  the  amount  that  xjj  can  be  increased  without  causing  an 
infeasibility  must  be  found.  This  is  the  subject  of  the  next  section. 


^  4.  The  Extended  Simplex  Method 

In  this  section  it  will  be  necessary  to  consider  situations  where  nonbasic  variables  are 
temporarily  positive,  as  in  the  Simplex  Method.  Assume  that  G  is  a  feasible  basis,  but  that 
the  corresponding  basic  feasible  solution  is  not  optimal  because  KTl  is  not  satisfied  for 

some  (I,J).  Then  the  reduced  cost  X]-p.jE(I,J)  is  negative,  so  xij  should  be  increased  as  in 
Manifold  suboptimization.  Since  the  generality  will  be  needed  later,  suppose  that  the 
nonbasic  variable  xjj  is  fixed  at  some  level  b  between  0  and  bj ,  and  that  the  balance 
equation  (6)  is  satisfied  for  each  component  T  of  G  by  z-1  when  xu=b.  If  the  nonbasic 

variable  xjj  is  increased  from  b  to  b+A,  let  the  solution  of  (6)  for  z  be  Z(T,A).  If  f()  is  a 
good  function,  Z(T,A)  is  a  continuous,  positive  function  of  A  with  Z(T,0)=1.  The 
corresponding  multipliers  (X(A),|i(A))  are  the  same  as  except  that  multipliers  for 

nodes  in  T  are  multiplied  by  Z(T,A).  Let  y(A)=g(|i(A)),  let  x(A)  be  the  unique  solution  of 
KT4-KT5  where  only  basic  variables  are  nonzero  except  that  nonbasic  xij=b-t-A.  let  A'  be 
the  largest  A  such  that  x(A)>(),  and  let  xk.l(A')=0.  In  other  words,  basic  variable  (K.L)  is 

the  one  first  driven  to  0  by  increasing  A.  If  x(A)>0  for  all  A^).  let  A'  and  (K.L)  be 
undefined.  If  (K.L)  is  defined  and  if  its  deletion  from  the  basis  breaks  any  cycle  that 
introducing  (I.J)  might  form,  then  the  next  basis  simply  replaces  (K,L)  with  (I.J).  (K.L)  is 
not  necessarily  defined  or  pan  of  a  cycle  even  if  it  is  defined,  but  Algorithm  A  below 
accounts  for  these  possibilities.  The  input  to  the  algorithm  is  a  feasible  basis  G^’, 
associated  multipliers  (>,^,)i^’),  and  a  nonbasic  edge  (1,J)  with  xij=()  and  negative  reduced 
cost.  The  output  is  a  different  feasible  basis  G'  at  least  as  good  as  G^’,  together  with 
associated  multipliers  (>,',|i’).  The  basic  idea  of  the  algorithm  is  to  increase  xu  in  steps 
until  finally  fhe  reduced  cost  is  zero.  If  I  and  .1  are  in  the  same  tree,  xu  is  increased  in  step 
A4  until  some  edge  is  deleted.  Eventually  this  operation  will  put  1  and  J  into  different  trees. 

When  I  and  J  are  in  different  trees,  either  an  edge  is  deleted  from  one  of  the  trees  (A'<A  in 
step  A.i)  or  else  (I.J)  is  finally  added  to  the  basis  in  the  temiinal  step  .Ad  Here  is 
Algorithm  A. 

Al)Letk=0,  XM) 

A2)  Let  G=Gk,  and  b=Xk. 

A3)  Let  Ti  be  the  tree  in  G^^  that  includes  I  and  Tj  be  the  tree  that  includes  J.  If  Tj^^Tj 
go  to  A5. 

A4)  Let  T=Ti.  Determine  A'  b>  solving  (6)  for  Z(T,A)  and  proceeding  as  in  the 
paragraph  above.  Let  (K.L)  be  the  basic  edge  in  T  whose  variable  is  first  driven  to  0,  let 
Gk^l=G>^-  (K.L),  (>,k+l,pk+l)=(X(^')  and  Xk+l=Xk+A'.  Go  to  A2. 

A5)  Let  T=T].  Determine  A'  by  solving  (6)  for  Z(T.A)  and  proceeding  as  in  the 
paragraph  above.  Let  (K.L)  be  the  basic  edge  in  T  whose  variable  is  first  driven  to  0,  Let 

*  (Ki,Li)=(K.L)  and  Ai=A'. 

A6)  Let  T=Tj.  Determine  A'  by  solving  (6)  for  Z(T,A)  and  proceeding  as  in  the 
►  paragraph  above.  If  (K.L)  is  undefined  let  Aj=oo;  otherwise,  let  (Kj,Lj)=(K,L)  and 
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Aj=A'. 

A”  t  I-et  A*  be  the  smallest  nonnegative  solution  of  the  equation 
X*iZ(  ri,A)=|i*jZ(Tj,A)E(I,J),  with  A*=oo  if  there  is  no  such  solution. 

A8)  Let  A’=min{Ai,Aj,A*}.  Let  Xi*^+*=A,i*^Z(Ti,A')  for  ieTj,  or  >.i'^'^*=>.i*^Z(Tj,A') 
for  ieTj,  or  otherwise  Define  ji*^'*’*  similarly.  Let  X^+'=X'^+A'. 

A9)  If  A’=A*,  set  G'=G'‘+(LJ)  and  Stop. 

AlO)  If  A'=Ai,  let  G''+l=Gk-  (Ki,Li)  and  k=k+l.  Go  to  A2. 

A 1 1)  If  A'=Aj,  let  1  ^G*'  -  (Kj.Lj)  and  k=k+ 1 .  Go  to  A2. 

Theorem  3;  Algorithm  A  stops  after  finitely  many  steps.  The  output  is  a  feasible  basis 
different  from  the  input,  with  associated  objective  function  at  least  as  good. 

Proof:  A’  and  (K,L)  are  well  defined  in  A4  and  A5  because  row  I  is  included  in  Tj  and 

therefore  A’  cannot  exceed  bi.  If  I  and  J  are  included  in  the  same  subtree,  then  the  reduced 

cost  is  negative  when  A=0  and  proportional  to  Z(T,A),  therefore  negative  for  all  A.  If  I 

and  J  are  in  different  trees,  then  the  reduced  cost  is  negative  for  A^*,  since  the  functions 

in  A7  are  continuous  in  A.  Since  A'^*  in  A8,  the  reduced  cost  is  negative  throughout  and 
there  is  no  possibility  that  the  objective  function  might  increase.  It  is  also  clear  that  G'  is 
different  from  G.  since  G'  includes  (I,J)  and  G  does  not.  TTie  only  question  is  whether  A9 
is  encountered  after  finitely  many  steps. 

NVlien  an  edge  is  deleted,  the  effect  on  the  basis  is  always  that  there  is  one  more  tree 
component  and  one  less  edge.  Since  A2  can  only  be  revisited  after  deleting  an  edge,  the 
number  of  edges  in  G*^  is  an  upper  bound  on  the  number  of  time.s  A2  is  encountered.  Since 
A2  is  encountered  as  frequently  as  any  other  step,  the  proof  is  complete.  QED 

The  pivoting  operation  described  above  differs  from  the  comparable  Linear 
Programming  operation  in  that  the  number  of  edges  that  have  to  leave  the  basis  in  order  to 
get  (l.J)  in  is  not  always  1;  it  may  be  any  nonnegative  integer.  It  differs  essentially  in  this 
respect  from  ZangwiU's  (1970)  Convex  Simplex  Method  (CSM).  In  the  CSM  the  number 
of  basic  variables  is  constant,  and  every  pivot  replaces  a  basic  variable  with  a  nonbasic 
variable.  When  a  local  minimum  is  encountered  in  CSM,  the  corresponding  variable 
remains  nonbasic,  but  positive,  and  consequently  there  is  no  unique  solution  associated 
with  a  given  basis.  In  the  method  discussed  above,  the  number  of  basic  variables  is  not 
constant,  but  the  idea  that  nonbasic  variables  are  zero  is  preserved,  along  with  the  one-to- 
one  relationship  between  bases  and  solutions.  TTiis  relationship  implies  that  the  number  of 
pivots  is  bounded  by  the  number  of  bases,  whereas  the  number  of  pivots  required  by  the 
CS.M  may  be  unbounded. 

5 .  Implementation  in  the  Exponential  Case 

Two  FORTRAN  programs  have  been  written  that  implement  the  methods  of  sections  3 
and  4  in  the  Exponential  case:  MANEFO  in  the  case  of  Manifold  suboptimization  and 
SIMPLX  for  the  Extended  Simplex  Method.  Each  of  these  programs  exploits  the  forest 
structure  of  the  basis  to  minimize  storage  and  facilitate  computations.  The  forest  is  actually 
represented  as  a  single  tree  by  introducing  an  "earth”  node  to  which  all  of  the  components 
of  the  basis  are  connected  through  a  root  node  in  each  tree;  however,  in  this  exposition  the 
term  "tree"  will  continue  to  mean  a  component  of  the  basis  as  heretofore.  Most  of  the 
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operations  associated  with  pivoting  are  similar  to  those  encountered  in  transshipment 
problems.  Since  these  are  well  discussed  in  Bradley,  Brown,  and  Graves  (1977),  only  the 
aspects  associated  with  the  nonlinear  nature  of  the  objective  function  will  be  discussed 
here. 

Consider  first  the  problem  of  deleting  an  edge  from  the  basis.  If  the  edge  is  in  tree  T, 
then  deleting  it  amounts  to  pruning  off  a  branch  from  T.  A  reduced  version  of  T  will 
remain,  along  with  a  new  tree  T'  consisting  of  the  pruned  off  branch.  Although  equation 
(1)  will  hold  in  all  trees  before  the  pruning,  it  generally  will  not  hold  in  either  T  or  T’ 

afterwards.  The  multipliers  (X,|i)  must  be  scaled  to  make  (1)  hold  in  T  and  T',  bearing  in 

mind  that  y  depends  on  |i  and  that  nodes  in  T  will  have  a  different  z-factor  than  nodes  in 
T'.  These  z-factors  can  be  computed  using  the  formula  given  in  Example  1  of  Section  2. 

The  computation  in  T  requires  the  sum  ^  fij.  This  number  could  be  computed  by 

jeC(T’) 

summing  |ij  over  the  appropriate  set,  but  in  fact  both  MANIFO  and  SIMPLX  utilize  an 

additional  array  SDM(k)  to  store  the  sum  of  |ij  over  the  set  consisting  of  k  and  all  of  its 
column  successors  in  whatever  (rooted)  tree  k  belongs  to.  Thus,  if  (I,J)  is  deleted  from  a 
tree  with  root  node  K,  and  if  (say)  J  is  the  predecessor  of  I,  then  the  required  sum  in  T'  is 
SDM(I)  and  the  required  sum  in  T  is  SDM(K)-SDM(I).  Calculating  the  z-factor  then 

requires  a  single  exponentiation  in  each  tree.  After  updating  (X,|i,y),  the  allocations  in  T 
and  T'  are  computed  as  described  in  part  b)  of  Theorem  2.;  a  preorder  traversal  array 
permits  this  to  be  done  in  one  pass.  The  allocations  themselves  are  stored  in  array  X,  with 
X(k)  representing  the  allocation  from  node  k  to  its  predecessor  if  k  is  a  row,  or  from  the 
predecessor  to  k  if  k  is  a  column.  All  arrays  are  of  length  n+m+l  in  both  programs  (node 
length);  in  fact,  the  only  function  requiring  mn  operations  is  that  of  determining  the  entering 

edge  by  searching  for  the  largest  ratio  )tjE(i,j)Ai . 

The  addition  of  edge  (I,J)  to  the  basis  requires  the  joining  of  one  tree  to  another, 
possibly  after  or  in  conjunction  with  another  pruning  operation.  Lf  T  containing  J  is  to  be 
joined  to  T.  T'  is  first  "rehung"  so  that  the  root  node  is  J.  Next  T  is  scaled  so  that 

Xi=ujE(I,J),  which  requires  a  logarithm  in  updating  the  potentials  in  T.  The  unified  tree  is 
then  scaled  as  in  the  previous  paragraph.  It  is  during  the  subsequent  calculations  of  X  that 
the  possibility  of  negative  allocations  must  be  provided  for  in  MANIFO. 

SIMPLX  uses  all  of  the  arrays  in  MANIFO  plus  one  additional.  The  additional  array 

W  is  associated  with  determining  how  big  A  can  be  before  some  basic  variable  becomes  0, 
as  will  be  explained  shortly.  In  the  Exponential  case  the  function  Z(T,A)  satisfies 
ln(Z(T,A))=D(T)A,  where 

D(T)=[Xi*5(T,I)-|aj*E(I,J)5(T,J)]/  ^  pj.  (7) 

j£C(T) 

Note  that  D(T)  is  easy  to  compute  because  the  denominator  is  already  stored,  and  also  that 
D(T)  will  be  negative  in  step  A5,  positive  in  A6,  and  negative  (since  the  reduced  cost  is 

negative)  in  A4.  Multiplying  Pj  by  Z(T,A)  is  equivalent  to  subtracting  ln(Z(T,A))  from  yj, 

so  yj(A)=yj(0)-D(T)A  for  jeC(T).  Since  this  is  a  linear  function  of  A,  allocation  xij(A)  will 

also  be  a  linear  function  of  A  for  all  edges  in  T.  The  array  element  W(k)  contains  the  rate  at 

which  the  allocation  associated  with  node  k  decreases  with  A.  The  calculation  of  W  is 
equivalent  in  difficulty  to  the  calculation  of  X  in  MANEFO.  Given  W,  it  is  simple  to  obtain 
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the  limiting  value  A'  by  comparing  XAV  ratios.  The  X  array  can  be  updated  once  A'  is 
known. 

Both  procedures  require  a  basis  for  a  starting  point.  This  is  obtained  by  first  including 
all  edges  of  the  form  (k,k).  If  there  are  more  columns  than  rows,  this  will  leave  some  trees 
consisting  of  singleton  columns.  If  there  are  more  rows  than  columns,  then  row  n+k  is 
assigned  to  column  k  for  l<k<n,  etc.,  until  finally  each  row  has  been  included  in  some 
edge.  The  resulting  basis  will  be  feasible,  but  unfortunately  some  edge  (ij)  for  which 
E(i,j)  =0  may  be  included.  The  simplest  (and  the  operative)  remedy  is  to  simply  change  * 

E(i  j)  to  some  small  positive  number,  since  in  any  case  either  procedure  will  quickly  delete 
such  an  edge. 

There  is  potential  for  roundoff  errors  to  accumulate  in  the  repeated  updating  of  floating 
point  arrays,  so  MANIFO  and  SIMPLX  both  use  double  precision  arithmetic.  Neither 
program  incorporates  any  protection  against  degeneracy,  but  so  far  there  have  been  no 
instances  of  cycling. 

To  give  the  reader  a  better  idea  of  the  actual  operation  of  the  algorithms.  Figure  1 
shows  the  sequence  of  bases  considered  by  MANIFO  in  the  process  of  solving  a  problem 
where  m=3  and  n=4.  The  basis  is  shown  as  a  single  tree  with  row  nodes  being 
distinguished  from  column  nodes  (as  they  are  in  MANIFO)  by  representing  row  nodes 
with  negative  numbers.  Node  0  is  "eanh".  With  the  convention  of  including  all  edges 
incident  to  node  0,  a  basis  always  has  exactly  n+m  edges.  Although  X  is  stored  as  a  node 
length  array  in  MANIFO,  the  allocations  are  shown  in  figure  1  in  the  more  familiar  matrix 
form.  There  are  (2^-l)^=3375  bases  for  this  problem,  of  which  MANIFO  considers  8 
feasible  and  2  nonfeasible  before  finding  one  that  is  optimal.  In  the  course  of  doing  this 
MANIFO  requires  18  logarithms  or  exponentials.  These  statistics  are  the  same  for 
SIMPLX  in  this  small  problem;  differences  emerge  only  in  problems  where  M.ANIFO 
introduces  more  than  1  negative  allocation  at  a  time. 

6 .  Computational  Comparisons 

MANIFO  and  SIMPLX  were  compared  on  randomly  generated  problems  to  determine 
which  program  is  fastest.  Problems  generated  were  such  that 
E(i,j)  is  uniform  [0,1] 
bj  is  uniform  [0,10] 

V(j)  is  uniform  [0,100] 

Eleven  (m,n)  pairs  were  tested;  (10,10),  (7,20),  (20,7),  (20,20),  (30,30),  40,40), 

(20,40),  (40,20),  (50,50),  (60,60),  and  (70,70).  Each  of  eleven  (E,b)  configurations  w'as 
combined  with  each  of  five  V  configurations  to  generate  a  problem.  For  each  of  these  55 
problems,  the  number  of  nonlinear  operations  (logarithms  and  exponentials)  NL  was 
recorded  along  with  the  CPU  time  T  in  seconds  required  to  set  up  the  initial  guess  and 
solve  the  problem,  m,  n,  NL,  and  T  were  then  transformed  by  taking  logarithms  and  a 
linetir  regression  performed. 

There  turns  out  to  be  very  little  difference  between  MANIFO  and  SIMPLX.  The  fitted 
formulas  are  in  each  case  approximately 

T  =  (.344  sec)(m/46)l -05(11/1 22)-8  exp(m/46+n/122)  (R2=99%)  (8) 

NL=1.4m^  (R2  =  98%)  (9) 

All  runs  were  made  on  the  Naval  Postgraduate  School's  IBM3033AP.  Since 
GAMS/MINOS  is  also  available  on  that  system,  in  one  instance  a  problem  with  m=20  and  t 

n=20  was  solved  using  GAMS/MINOS.  The  time  required  was  1.88  seconds  ("MINOS 
TIME"  in  the  output),  which  is  47  times  as  large  as  the  time  required  by  MANIFO  (.04 
set.ond.s)  for  the  same  problem.  Further  comparisons  were  not  made  because  MINOS  is  a  , 

general  purpose  solver  that  does  not  exploit  the  network  nature  of  the  constraints.  A  bener 
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comparison  would  be  with  the  GENOS  network  optimizer  (Mulvey  and  Zlemos  (1987)). 
Attempts  to  accomplisii  this  are  in  progress,  but  so  far  the  lack  of  an  exponential  function  in 
GENOS  1 .0  has  proved  to  be  a  roadblock. 

The  reason  for  keeping  track  of  the  number  of  exponentials  and  logarithms  NL  in 
running  MANIFO/SINIPLX  is  that  it  was  initially  anticipated  that  such  nonlinear  operations 
would  require  relatively  large  amounts  of  time.  This  turns  out  not  to  be  the  case.  When 
(m,n)  =  (70,70),  formulas  (8)  and  (9)  predict  T  =  .28  seconds  and  NL  =  820.  Since  an 
IBM3033AP  requires  only  about  4p.sec  for  a  double  precision  logarithm  or  exponential,  the 
amount  of  time  spent  in  nonlinear  operations  is  only  about  (820)(4xl0‘^)  =  .0033  secords; 
MANEFO  and  SIMPLX  are  each  mainly  occupied  with  manipulating  and  scaling  the 
various  arrays  required  to  represent  the  current  basic  solution  as  a  tree,  rather  than  solving 
the  balance  equation.  This  is  encouraging  for  applications  where  the  objective  function, 
although  "goad,"  might  require  a  numeric  rather  than  analytic  solution  of  the  balance 
equation. 
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Figure  1 

Solution  of  a  3x4  Problem  Using  MANIFO  (Continued) 
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